Fig.1 - Date of 1st lab results

Weekly temporal distribution of first positive COVID-19 lab results throughout the study period.

FirstDate <- 
  BD_chusj %>% 
  select(DateOfFirstPositiveLabResult) %>% 
  group_by(date = DateOfFirstPositiveLabResult) %>% 
  summarise(freq = n())

FirstDate_daily   <- xts::as.xts(FirstDate$freq, order.by = as.Date(FirstDate$date))
FirstDate_weekly  <- xts::apply.weekly(FirstDate_daily, sum)

xts::plot.xts(FirstDate_weekly, col = "cadetblue4",
              lwd = 3, main = "Date of first positive lab result - weekly",
              yaxis.right = FALSE)








Fig.2 - Radar Plots

Radar plots representing medication usage in subgroups: (A) ICU-admitted patients, (B) type of hospital discharge, and (C) type of MRCS. In each radar, the lines represent distinct patient categories, and the distance from each point to the origin indicates the percentage of patients using drugs from a specific ATC group. For instance, a line extending to the 75% marker under the label “ATC.A” suggests that 75% of those patients are being medicated with drugs from the ATC.A group. It’s important to note that comparisons between lines within the same radar provide insights into discrepancies in medication usage among the subgroups; however, specific values and detailed interpretation should be derived from the study’s context.

#' Define function to revert order of variables in radar plot
clockwise_order <- function(ds){
  # ds <-ATCs.lv1
  ds_rev <- rev(ds)
  ds_rev <- ds_rev[, c(ncol(ds_rev), 1:(ncol(ATCs.lv1)-1))]
  return(ds_rev)
}


# Define multiframe
par(mfrow = c(1, 3))
legend.x = -1.3
legend.y = -0.9


# Fig.2A - ATCs VS Intensive Care -----------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(IntensiveCare = BD_chusj$IntensiveCare,
                  ATCs.lv1)

IC_yes_count <- ATCs.lv1$IntensiveCare[ATCs.lv1$IntensiveCare == 1] %>% length()
IC_no_count  <- ATCs.lv1$IntensiveCare[ATCs.lv1$IntensiveCare == 0] %>% length()

# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
  # Intensive Care == YES
  IC_yes = ATCs.lv1 %>%
    filter(IntensiveCare == 1) %>%
    {sapply(., sum) / IC_yes_count} %>%
    round(digits = 3) %>% t() %>% as.data.frame(),
  # Intensive Care == NO
  IC_no = ATCs.lv1 %>%
    filter(IntensiveCare == 0) %>%
    {sapply(., sum) / IC_no_count} %>%
    round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$IntensiveCare <- NULL

# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)

# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
                  min = rep(0, length(ATCs.lv1)),
                  ATCs.lv1)

# Define colors
colors <- c("indianred2", "forestgreen")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.2)

# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)

radarchart(ATCs.lv1,
           axistype = 1,
           # custom polygon
           pcol = colors_border,  # line color
           pfcol = colors_in,     # fill-in color
           plwd = 2,              # line width
           plty = 1,              # line type
           # custom the grid
           cglcol = "grey20",     # net color
           cglty = 3,             # net line type
           axislabcol = "grey50", # net label color
           cglwd = 0.8,           # net width
           caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
           # custom labels
           calcex = 0.8,
           vlcex = 0.7,
           title = "Relative use of medication in ICU-admitted\nand non-ICU-admitted patients",
           sub = paste0("ICU-admitted patients: n=", IC_yes_count, "\n",
                        "Non-ICU-admitted patients: n=", IC_no_count)
)

legend(x = legend.x,
       y = legend.y,
       legend = rownames(ATCs.lv1[-c(1,2),]),
       bty = "n",
       pch = 20,
       col = colors_border,
       text.col = "grey40",
       cex = 0.9,
       pt.cex = 2)

text(-1.2, 1.2, "(A)", cex = 1.4, pos = 1, col = "black")



# Fig.2B - ATCs vs Hospital Discharge -------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(Outcome = BD_chusj$Outcome,
                  ATCs.lv1)

Outcome_Recovered_count    <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 0] %>% na.omit() %>% length()
Outcome_Transferred_count  <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 2] %>% na.omit() %>% length()
Outcome_Deceased_count     <- ATCs.lv1$Outcome[ATCs.lv1$Outcome == 1] %>% na.omit() %>% length()

# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
  Recovered    = ATCs.lv1 %>% filter(Outcome == 0) %>% {sapply(., sum) / Outcome_Recovered_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
  Transferred  = ATCs.lv1 %>% filter(Outcome == 2) %>% {sapply(., sum) / Outcome_Transferred_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
  Deceased     = ATCs.lv1 %>% filter(Outcome == 1) %>% {sapply(., sum) / Outcome_Deceased_count} %>% round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$Outcome <- NULL

# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)

# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
                  min = rep(0, length(ATCs.lv1)),
                  ATCs.lv1)

# Define colors
colors <- c("forestgreen", "dodgerblue3", "firebrick3")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.1)

# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)

radarchart(ATCs.lv1,
           axistype = 1,
           # custom polygon
           pcol = colors_border,  # line color
           pfcol = colors_in,     # fill-in color
           plwd = 2,              # line width
           plty = 1,              # line type
           # custom the grid
           cglcol = "grey20",   # net color
           cglty = 3,  # net line type
           axislabcol = "grey50", # net label color
           cglwd = 0.8,         # net width
           caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
           # custom labels
           calcex = 0.8,
           vlcex = 0.7,
           title = "Relative use of medication in patients\n based on hospital discharge",
           sub = paste0("Patients Recovered: n=", Outcome_Recovered_count, "\n",
                       "Patients Transfered: n=", Outcome_Transferred_count, "\n",
                       "Patients Deceased: n=", Outcome_Deceased_count)
)

legend(x = legend.x,
       y = legend.y,
       legend = rownames(ATCs.lv1[-c(1,2),]),
       bty = "n",
       pch = 20,
       col = colors_border,
       text.col = "grey40",
       cex = 0.9,
       pt.cex = 2)

text(-1.2, 1.2, "(B)", cex = 1.4, pos = 1, col = "black") 





# Fig.2C - ATCs vs Respiratory Support ------------------------------------
ATCs.lv1 <- BD_chusj[, c(469, 378:390)]
ATCs.lv1 <- cbind(BD_chusj[25:28],
                  ATCs.lv1)

Support_ECMO_count <- ATCs.lv1$ECMO[ATCs.lv1$ECMO == 1] %>% na.omit() %>% length()
Support_HFO_count <- ATCs.lv1$HFO[ATCs.lv1$HFO == 1] %>% na.omit() %>% length()
Support_NIV_count <- ATCs.lv1$NIV[ATCs.lv1$NIV == 1] %>% na.omit() %>% length()
Support_IMV_count <- ATCs.lv1$IMV[ATCs.lv1$IMV == 1] %>% na.omit() %>% length()

# Calculate the percentage of patients in and out of IC for each ATC group
ATCs.lv1 <- rbind(
  ECMO = ATCs.lv1 %>% filter(ECMO == 1) %>% {sapply(., sum) / Support_ECMO_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
  HFO = ATCs.lv1 %>% filter(HFO == 1) %>% {sapply(., sum) / Support_HFO_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
  NIV = ATCs.lv1 %>% filter(NIV == 1) %>% {sapply(., sum) / Support_NIV_count} %>% round(digits = 3) %>% t() %>% as.data.frame(),
  IMV = ATCs.lv1 %>% filter(IMV == 1) %>% {sapply(., sum) / Support_IMV_count} %>% round(digits = 3) %>% t() %>% as.data.frame()
)
ATCs.lv1$ECMO <- ATCs.lv1$HFO <- ATCs.lv1$NIV <- ATCs.lv1$IMV <- NULL

# Find max and round it (ceiling) to the nearest .1 multiple
max.axis <- plyr::round_any(max(ATCs.lv1), 0.1, f = ceiling)

# Prepare data for radar chart
ATCs.lv1 <- rbind(max = rep(max.axis, length(ATCs.lv1)),
                  min = rep(0, length(ATCs.lv1)),
                  ATCs.lv1)

# Define colors
colors <- c("firebrick2", "dodgerblue2", "forestgreen", "darkorange")
colors_border <- colors
colors_in <- scales::alpha(colors, alpha = 0.1)

# Revert order of columns
ATCs.lv1 <- clockwise_order(ATCs.lv1)

radarchart(ATCs.lv1,
           axistype = 1,
           # custom polygon
           pcol = colors_border,  # line color
           pfcol = colors_in,     # fill-in color
           plwd = 2,              # line width
           plty = 1,              # line type
           # custom the grid
           cglcol = "grey20",   # net color
           cglty = 3,  # net line type
           axislabcol = "grey50", # net label color
           cglwd = 0.8,         # net width
           caxislabels = scales::percent(seq(0, max.axis, length.out = 5)),
           # custom labels
           calcex = 0.8,
           vlcex = 0.7,
           title = "Relative use of medication in patients based on\n Mechanical Respiratory and Circulatory Support",
           sub = paste0("Patients using support:\n",
                        "ECMO: n=", Support_ECMO_count, ";  ",
                        "HFO: n=", Support_HFO_count, ";\n",
                        "NIV: n=", Support_NIV_count, ";  ",
                        "IMV: n=", Support_IMV_count
           )
)

legend(x = legend.x,
       y = legend.y,
       legend = rownames(ATCs.lv1[-c(1,2),]),
       bty = "n",
       pch = 20,
       col = colors_border,
       text.col = "grey40",
       cex = 0.9,
       pt.cex = 2)

text(-1.2, 1.2, "(C)", cex = 1.4, pos = 1, col = "black")

# Clear multiframe
par(mfrow = c(1,1))








Fig.3 - ATC families, 1st level

Heatmap illustrating co-frequencies of exposure among 1st level ATC drug groups. Each cell represents the relative co-frequency of patients using two drug groups concurrently. For instance, the cell “ATC.A x ATC.C”, with a value of 0.360, indicates that 36% of the total patients (n=2688) are simultaneously prescribed drugs from ATC groups A and C. The diagonal line, like the cell “ATC.C x ATC.C” showing 0.412, can be interpreted as the relative frequency of the use of that specific ATC drug group on its own, signifying that 41.2% of the total patients are prescribed only drugs from ATC group C. It’s important to note that the heatmap is mirrored, meaning data below the diagonal mirror those above it. Therefore, the cell “ATC.A x ATC.C” has the same value as the cell “ATC.C x ATC.A”.

#' 469 + 378-390: ATC codes - first level
ATC_lv1 <- BD_chusj[, c(469, 378:390)]

# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)


for (i in 1:length(ATC_lv1)) {
  for (j in i:length(ATC_lv1)) {
    heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
  }
}; rm(i, j)


heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value") 

# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)

# Plot heapmap
ggplot(heatmap.data, aes(x, y, fill = value)) + 
  geom_tile(color = "grey95",
            lwd = 1,
            linetype = 1) + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
          subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
  geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
  # scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
  scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
  guides(fill = guide_colourbar(barwidth = 0.5,
                                barheight = 20,
                                title = "%"))



Triangular version 1

Excluding diagonal line

# Reshape data to triangular matrix - version_1
heatmap.data_aux <- reshape2::dcast(heatmap.data, x ~ y)
n <- nrow(heatmap.data_aux)
mat <- matrix(NA, n, n)
mat[row(mat) + col(mat) <= n] <- 1
mat <- mat |> as.data.frame() |> rev()
data_vector <- as.vector(as.matrix(mat)) * as.vector(as.matrix(heatmap.data_aux[, 2:(n+1)]))
heatmap.data_diagonal_v1 <- cbind(heatmap.data_aux$x, as.data.frame(matrix(data_vector, nrow = n)))
names(heatmap.data_diagonal_v1) <- names(heatmap.data_aux)
heatmap.data_diagonal_v1 <- reshape2::melt(heatmap.data_diagonal_v1, value.name = "value", id.vars = "x")
names(heatmap.data_diagonal_v1) <- c("x", "y", "value")
heatmap.data_diagonal_v1 <- heatmap.data_diagonal_v1[complete.cases(heatmap.data_diagonal_v1),]


# Plot triangular heapmap v1
ggplot(heatmap.data_diagonal_v1, aes(x, y, fill = value)) + 
  theme_light() +
  geom_tile(color = "grey95",
            lwd = 1,
            linetype = 1) + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_line(color='white')) +
  ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
          subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
  geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
  # scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
  scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
  guides(fill = guide_colourbar(barwidth = 0.5,
                                barheight = 20,
                                title = "%"))



Triangular version 2

Excluding diagonal line

# Reshape data to triangular matrix - version_2
heatmap.data_aux <- reshape2::dcast(heatmap.data, x ~ y)
n <- nrow(heatmap.data_aux)
mat <- matrix(NA, n, n)
mat[row(mat) + col(mat) > n+1] <- 1
mat <- mat |> as.data.frame() |> rev()
data_vector <- as.vector(as.matrix(mat)) * as.vector(as.matrix(heatmap.data_aux[, 2:(n+1)]))
heatmap.data_diagonal_v2 <- cbind(heatmap.data_aux$x, as.data.frame(matrix(data_vector, nrow = n)))
names(heatmap.data_diagonal_v2) <- names(heatmap.data_aux)
heatmap.data_diagonal_v2 <- reshape2::melt(heatmap.data_diagonal_v2, value.name = "value", id.vars = "x")
names(heatmap.data_diagonal_v2) <- c("x", "y", "value")
heatmap.data_diagonal_v2 <- heatmap.data_diagonal_v2[complete.cases(heatmap.data_diagonal_v2),]


# Plot triangular heapmap v2
ggplot(heatmap.data_diagonal_v2, aes(x, y, fill = value)) + 
  theme_light() +
  geom_tile(color = "grey95",
            lwd = 1,
            linetype = 1) + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_line(color='white')) +
  ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
          subtitle = paste0("Expressed as % of total patients (n = ", nrow(ATC_lv1)," patients)")) +
  geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
  # scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
  scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
  guides(fill = guide_colourbar(barwidth = 0.5,
                                barheight = 20,
                                title = "%"))








Fig.4 - ATC groups 2nd level

Heatmaps illustrating co-frequencies among 2nd level ATC drug groups for the most prevalent associations from the 1st level ATC drug group heatmap: (A) ATC A and B, (B) ATC A and N, (C) ATC B and N, and (D) ATC B and J. Each cell represents the relative co-frequency of patients using two drug groups concurrently. For instance, in the heatmap (A), the cell “ATC.B01 x ATC.A02” with a value of 0.333 indicates that 33.3% of the total patients (n=2688) are concurrently prescribed drugs from 2nd level ATC groups B01 and A02.

heatmap_2groups <- function(index1 = NA,
                            index2 = NA,
                            verbose = TRUE) {
  
  #' This function plots a heatmap of 2 ATC groups at 2nd level
  #' inputs: the indexes on the 2nd levels codes for each group
  #'   ATC.A: 362:377 (16 codes)
  #'   ATC.B: 391:395 (5 codes)
  #'   ATC.J: 425:430 (6 codes)
  #'   ATC.N: 441:447 (7 codes)
  
  ATC_group1 <- BD_chusj[, index1]
  ATC_group2 <- BD_chusj[, index2]
  
  # Create an empty dataset
  heatmap.matrix <- data.frame(matrix(NA,
                                      ncol = ncol(ATC_group1),
                                      nrow = ncol(ATC_group2)))
  colnames(heatmap.matrix) <- names(ATC_group1)
  rownames(heatmap.matrix) <- names(ATC_group2)
  
  ATC_groups12 <- cbind(ATC_group1, ATC_group2)
  for (g1 in seq_along(ATC_group1)) {
    for (g2 in seq_along(ATC_group2)) {
      heatmap.matrix[g2, g1] <- 
        nrow(ATC_groups12[, c(g1, g2 + 
                                ncol(ATC_group1))][ATC_groups12[,g1] == 1 &
                                                     ATC_groups12[,g2 + ncol(ATC_group1)] == 1,])
    }
  }
  
  heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
  colnames(heatmap.data) <- c("x", "y", "value")
  
  # change values to percentage
  heatmap.data$value <- round(heatmap.data$value / nrow(BD_chusj), digits = 3)
  
  family.group1 <- substring(names(ATC_group1)[1], 5, 5)
  family.group2 <- substring(names(ATC_group2)[1], 5, 5)
  
  p <- 
    ggplot(heatmap.data, aes(x, y, fill = value)) + 
    geom_tile(color = "grey95",
              lwd = 1,
              linetype = 1) +
    theme(axis.text.x = element_text(angle = 0),
          axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    ggtitle(label = paste("Relative co-frequences in 2nd-level ATC drug groups",
                          family.group1, "and", family.group2),
            subtitle = "Expressed as % of total patients (n = 2688)") +
    geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
    # scale_fill_gradient(low = "lightsteelblue1", high = "darkolivegreen3") +
    scale_fill_gradient(low = "#F9FFFD", high = "#0C8A00") +
    # coord_fixed() +
    guides(fill = guide_colourbar(barwidth = 0.5,
                                  barheight = 7,
                                  title = "%"))

  # Find most prevalent combinations
  if (verbose) {
    heatmap.data %>%
      arrange(desc(value)) %>%
      mutate(perc = round(value / nrow(BD_chusj), 3)) %>%
      relocate(y) %>%
      rename(!!family.group1 := y, !!family.group2 := x, freq = value) %>%
      head() %>% print()
  }
  
  return(p)
}
# Plot all heatmaps together
p1 <- heatmap_2groups(index1 = 362:377, index2 = 391:395, verbose = F) # A.B
p2 <- heatmap_2groups(index1 = 391:395, index2 = 441:447, verbose = F) # B.N
p3 <- heatmap_2groups(index1 = 391:395, index2 = 425:430, verbose = F) # B.J
p4 <- heatmap_2groups(index1 = 362:377, index2 = 441:447, verbose = F) # A.N

cowplot::plot_grid(p1, p4, p2, p3,
                   ncol = 2, 
                   nrow = 2,
                   rel_heights = c(2,1),
                   scale = 0.9,
                   labels = c("(A)", "(B)", "(C)", "(D)"),
                   label_size = 14
)








Fig.5 - ICU-admittance & Deceased

Heatmaps illustrating co-frequencies of exposure among 1st level ATC drug groups for (A) ICU-admitted patients and (B) deceased patients. Each cell represents the relative co-frequency of patients using two drug groups concurrently. Taking the heatmap (A) as an example, the cell “ATC.B x ATC.A”, with a value of 0.934, indicates that 93.4% of ICU-admitted patients (n=858) are simultaneously prescribed drugs from ATC groups B and A. The diagonal cells, such as the cell “ATC.B x ATC.B” showing 0.997, can be interpreted as the relative frequency of the use of that specific ATC drug group alone, signifying that 99.7% of ICU-admitted patients are prescribed only drugs from ATC group B. It’s important to note that the heatmap is mirrored, meaning data below the diagonal mirror those above it. Therefore, the cell “ATC.A x ATC.B” has the same value as the cell “ATC.B x ATC.A”.

# Fig.5A ------------------------------------------------------------------
ATC_lv1 <- BD_chusj[, c(469, 378:390)]
ATC_lv1 <- cbind(IntensiveCare = BD_chusj$IntensiveCare,
                 ATC_lv1)

# Filter only those patients in Intensive Care
ATC_lv1 <- ATC_lv1 %>% filter(IntensiveCare == 1)
ATC_lv1$IntensiveCare <- NULL

# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)

for (i in 1:length(ATC_lv1)) {
  for (j in i:length(ATC_lv1)) {
    heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
  }
}; rm(i, j)

heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value") 

# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)

# store graphic
fig.5A <-
  ggplot(heatmap.data, aes(x, y, fill = value)) + 
  geom_tile(color = "grey95",
            lwd = 1,
            linetype = 1) + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
          subtitle = paste0("Expressed as % of patients in Intensive Care (n = ", nrow(ATC_lv1), " patients)")) +
  geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
  scale_fill_gradient(low = "#F9FFFD", high = "#a66e29") + #8068b0 #a66e29
  guides(fill = guide_colourbar(barwidth = 0.5,
                                barheight = 20,
                                title = "%"))





# Fig.5B ------------------------------------------------------------------
ATC_lv1 <- BD_chusj[, c(469, 378:390)]
ATC_lv1 <- cbind(Outcome = BD_chusj$Outcome,
                 ATC_lv1)

# Filter only those patients in Intensive Care
ATC_lv1 <- ATC_lv1 %>% filter(Outcome == 1)
ATC_lv1$Outcome <- NULL

# Create an empty dataset
heatmap.matrix <- data.frame(matrix(NA, ncol = length(ATC_lv1), nrow = length(ATC_lv1)))
colnames(heatmap.matrix) <- rownames(heatmap.matrix) <- names(ATC_lv1)

for (i in 1:length(ATC_lv1)) {
  for (j in i:length(ATC_lv1)) {
    heatmap.matrix[i,j] <- heatmap.matrix[j,i] <- nrow(ATC_lv1[,c(i,j)][ATC_lv1[,i] == 1 & ATC_lv1[,j] == 1,])
  }
}; rm(i, j)

heatmap.data <- reshape2::melt(as.matrix(heatmap.matrix))
colnames(heatmap.data) <- c("x", "y", "value") 

# change values to percentage
heatmap.data$value <- round(heatmap.data$value / nrow(ATC_lv1), digits = 3)

# store graphic
fig.5B <-
  ggplot(heatmap.data, aes(x, y, fill = value)) + 
  geom_tile(color = "grey95",
            lwd = 1,
            linetype = 1) + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  ggtitle(label = "Relative co-frequencies in 1st-level ATC drug groups",
          subtitle = paste0("Expressed as % of deceased patients (n = ", nrow(ATC_lv1), " patients)")) +
  geom_text(aes(label = sprintf("%0.3f", value)), color = "gray30", size = 3) +
  scale_fill_gradient(low = "#F9FFFD", high = "#cc3746") +
  guides(fill = guide_colourbar(barwidth = 0.5,
                                barheight = 20,
                                title = "%"))



# Plot all heatmaps together
cowplot::plot_grid(fig.5A, fig.5B,
                   ncol = 2, 
                   nrow = 1,
                   # rel_heights = c(2,1),
                   scale = 0.9,
                   labels = c("(A)", "(B)"),
                   label_size = 14
)